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1. Introduction 

The last decade has witnessed a tremendous amount of research efforts in the physics of 
atomic Bose- Einstein condensates (BECs) [H [2]. The study of BECs has yielded a wide 
array of interesting phenomena, not only because of very precise experimental control 
that exists over the relevant experimental procedures [3j H], but also because of the 
intimate connections of the description of dilute-gas BECs with other areas of physics, 
such as superfluidity, superconductivity, lasers and coherent optics, nonlinear optics, 
and nonlinear wave theory. Of particular emphasis in much experimental and theoretical 
work is the setting of a BEC trapped in periodic potentials, usually combined with an 
additional harmonic trapping potential. From the standpoint of nonlinear interactions, 
mathematical descriptions of BECs held in purely harmonic traps are now well known. 
Nevertheless, apart from studies focusing on the transition between superfluidity and an 
insulator state [5], relatively little attention has been given to an understanding of the 
varieties of many-body states that may possibly exist with intermediate lattice strengths, 
where phase coherence is maintained across the sample. A more complete understanding 
of BEC behavior in such lattice potentials is relevant and important to current work with 
BECs, and to an even broader array of topics, in particular discrete nonlinear optics 
and nonlinear wave theories. Such regimes of BEC physics are experimentally and 
theoretically accessible, and comparisons between theoretical and experimental results 
are certainly possible. Here, we present a theoretical stability examination of BECs 
with either attractive or repulsive interatomic interactions in a combined harmonic and 
periodic potential. 

Many of the common elements between BECs and other areas of physics, and in 
particular optics, originate in the existence of macroscopic coherence in the many-body 
state of the system. Mathematically, BEC dynamics are therefore often accurately 
described by a mean-field model, namely a partial differential equation of the nonlinear 
Schrodinger (NLS) type, the so-called Gross-Pitaevskii (GP) equation [TJ [2], [3]. The 
GP equation is particularly successful in drawing connections between BEC physics, 
nonlinear optics and nonlinear wave theories, with vortices and solitary waves examples 
of common elements between these areas. The GP equation is a classical nonlinear 
evolution equation (with the nonlinearity originating from the interatomic interactions) 
and, as such, it permits the study of a variety of interesting nonlinear phenomena. 
These phenomena have primarily been studied by treating the condensate as a purely 
nonlinear coherent matter- wave, i.e., from the viewpoint of the nonlinear dynamics of 
solitary waves. Relevant studies have already been summarized in various books (see, 
e.g., Ref. [I]) and reviews (see, e.g., Refs. [6] for bright matter-wave solitons, [HE] for 
vortices in BECs, [9] for dynamical instabilities in BECs, [TQl [11] for nonlinear dynamics 
of BECs in optical lattices). 

On the other hand, many static and dynamic properties of BECs confined in various 
types of external potentials can be studied by starting from the non- interacting limit, 
where the nonlinearity is considered to be negligible. The basic idea of such an approach 
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is that in the absence of interactions the GP equation is reduced to a linear Schrodinger 
equation for a confined single-particle state; in this limit, and in the case of, e.g., a 
harmonic external potential, the linear problem becomes the equation for the quantum 
harmonic oscillator characterized by discrete energies and corresponding eigenmodes 
|12j . Exploiting this simple physical picture, one may then use analytical and/or 
numerical techniques for the continuation of these linear eigenmodes supported by the 
particular type of the external trapping potential into nonlinear states as the interactions 
become stronger. This idea has been explored at the level of one- dimensional (ID) [13] 
and higher- dimensional states [14J in the case of a harmonic trapping potential, where 
nonlinear stationary modes were found from a continuation of the (linear) states of the 
quantum harmonic oscillator. The same problem has been studied in the framework 
of the so-called Feshbach resonance management technique in Ref. [15], where a linear 
temporal variation of the nonlinearity was considered. The continuation of the linear 
states to their nonlinear counterparts has also been explored from the point of view 
of bifurcation and stability theory [16]. Finally, in the same spirit but in the two- 
dimensional (2D) setting, radially symmetric nonlinear states of harmonically trapped 
"pancake-shaped" condensates were recently investigated in Ref. [17] . 

Importantly, all of the above studies provide a clear physical picture of how 
genuinely nonlinear states of harmonically confined BECs (such as dark and bright 
matter- wave solitons in ID or ring solitons and vortices in 2D), are connected to 
and emanate from the eigenmodes of the quantum harmonic oscillator. Similar 
considerations also hold for BECs confined in optical lattices. In this case, pertinent 
nonlinear stationary states (such as spatially extended nonlinear Bloch waves, truncated 
nonlinear Bloch waves, matter- wave gap solitons in ID, and gap vortices in 2D and 
3D), can be understood by the structure of the band-gap spectrum of the linear Bloch 
waves supported in the non-interacting limit (see, e.g., Ref. [18] and references therein). 
However, there are only few studies for condensates confined in both harmonic and 
optical lattice potentials, and these are basically devoted to the dynamics of particular 
nonlinear structures (such as dark [19J and bright [20] solitons in ID, and vortices in 2D 
[2T]). Thus, the structure of condensates confined in such superpositions of harmonic 
and periodic potentials remains, to the best of our knowledge, largely unexplored. 

Nevertheless, such a study is particularly relevant to current work with BECs, and 
even suggests new avenues for exploration. In particular, one might ask whether the 
addition of a weak optical lattice might increase the stability of excited states that 
are known to be unstable in harmonic traps. Stability, if it is found, may add new 
realistic options for the engineering of new quantum states of BECs. It may also 
be interesting to investigate the Mott insulator transition by starting from a stable 
excited state of a weak optical lattice. Also, the transport of excited states (which 
is not discussed in this paper) through a lattice structure may have application in 
future BEC interferometry experiments. Finally, the advances of far-off-resonant optical 
trapping techniques allow for the creation of strongly pancake-shaped condensates that 
may be confined by harmonic and spatially periodic components, and we expect that 
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the theoretical considerations described here may be directly explored using current 
experimental techniques. 

Our aim in the present work is to contribute to this direction and study the structure 
and the stability of a pancake-shaped condensate confined by the combination of a 
harmonic trap and a periodic potential, with periodicity in two orthogonal directions. 
We will adopt the above mentioned approach of continuation of linear states to nonlinear 
ones, thus providing a host of interesting solutions that have not been explored 
previously and yet should be tractable within presently available experimental settings. 
In particular, our analysis starts by first considering the non-interacting limit. In this 
regime, we employ a multiscale perturbation method (which uses the harmonic trap 
strength as a formal small parameter) to find the discrete energies and the corresponding 
eigenmodes of the pertinent single-particle Schrodinger equation with the combined 
harmonic and periodic potential. We then use this linear limit as a starting point 
for initializing a 2D "nonlinear solver" that identifies the relevant stationary nonlinear 
eigenstates as a function of the chemical potential (i.e., the nonlinearity strength) and 
of the optical lattice depth. 

Once the basic structure of the condensate is found, we subsequently perform 
a linear stability analysis of the nonlinear modes that can be initiated by the non- 
interacting ground and first few excited states. When nonlinear states are found to 
be unstable, we use direct numerical simulations to study their dynamics and monitor 
the evolution of the relevant instability. Essential results that will be presented below 
are the following: for a fixed harmonic trap strength, there exist certain regions in the 
parameter plane defined by the chemical potential and the optical lattice depth, where 
not only the ground state, but also excited states are stable or only weakly unstable. 
Particularly, an excited state with a shape resembling an out-of-phase matter-wave 
soliton pair (for attractive interactions) is found to persist for long times, being stable 
(weakly unstable) for attractive (repulsive) interatomic interactions. Thus, the ground 
state and the aforementioned excited state have a good chance to be observed in a 
real experiment with either attractive or repulsive pancake BECs. Similar conclusions 
can be drawn for more complex states, such as a quadrupolar one which may also be 
stable in the attractive case, however, higher excited states are typically more prone to 
instabilities, as is shown in our detailed numerics below. 

The paper is organized as follows: In Section II we present the model and study 
analytically the non-interacting regime. The continuation of the linear states to the 
nonlinear ones, as well as the stability properties of the nonlinear states are presented 
in Section III. Finally, in Section IV, we summarize our findings and present our 
conclusions. 

2. The model and its analytical consideration 

At sufficiently low temperatures, and in the framework of the mean-field approach, 
the condensate dynamics can be described by the order parameter \I/(r, t). We assume 
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that the condensate is kept in a highly anisotropic trap, with the transverse (x, y) and 
longitudinal (z) trapping frequencies chosen so that u x = uj y = u±_ <C u z . In such a case, 
the condensate has a nearly planar, so-called "pancake" shape (see, e.g., Refs. [22] for 
relevant experimental realizations), which allows us to assume a separable wave function, 
^ = Q(z)ip(x, y), where &{z) is the ground state of the respective quantum harmonic 
oscillator. Then, averaging of the underlying three-dimensional (3D) GP equation along 
the longitudinal direction z [23] leads to the following 2D GP equation for the transverse 
component of the wave function (see also Refs. [HI El H]): 

h 2 

ikdti) = -— + 92d\^\ 2 ^ + V ext (x, y)$. (1) 

Here, A = d 2 + d 2 is the 2D Laplacian, m is the atomic mass, and g2D — 9zdI '\/2Tia z is 
an effective 2D coupling constant, where g^D = Aith 2 a/m (a being the scattering length) 
and a z = yh/ muj z is the longitudinal harmonic oscillator length. Finally, the potential 
V ext (x, y) in the GP Eq. (JTJ is assumed to consist of a harmonic component and a square 
2D optical lattice (OL) created by two pairs of interfering laser beams of wavelength A: 

Vtxtfay) = ^mu 2 j_r 2 + V [cos 2 (kx) + cos 2 (ky)}. 

= V H (r) + V OL (x,y) (2) 

In the above expression, r 2 = x 2 + y 2 , while the optical lattice is characterized by two 
parameters, namely its depth Vq and its periodicity d = n/k = (X/2)/ sin(#/2), where 9 
is the angle between the two beams that create the x-direction lattice, and between the 
two beams that create the ^/-direction lattice. 

Measuring length in units of = d/ir, time in units of ujj} = H/El, and energy in 
units of El = 2E Tec = h 2 jma 2 L (where E Tec is the lattice recoil energy), the GP Eq. (TjQ) 
can be put into the following dimensionless form: 

idt1> = + s|VfV> + V(x, y)tp. (3) 

In the normalized GP Eq. ([3]), the wavefunction is rescaled as ip — > 
\J\g2D I / Eli[) exp [i(Vo / Ei)t] , the parameter s is given by s = sign(g 2 D) = ±1 (with 
s = +1 or s = - 1 corresponding, respectively, to repulsive or attractive interatomic 
interactions), while the normalized trapping potential V(x,y) is now given by: 

V{x, y) = ifiV + 1/ (cos(2x) + cos(2j/)). (4) 

In the above equation, the normalized lattice depth Vq is measured in units of 4£' rec , 
while the normalized harmonic trap strength is given by 

£1=4 = ^. (5) 

a± uj l 



where a± = yh/muj± is the transverse harmonic oscillator length. Note that using 
realistic parameter values (see, e.g., Ref. [24J), namely, a lattice periodicity 0.3 /im, a 
recoil energy E rec /h ~ 6KHz (assuming an atomic mass corresponding to 87 Rb), and 
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Figure 1. (Color Online) Panel (a) shows the energy spectra corresponding to 
a purely parabolic potential (pluses), a parabolic and lattice one found numerically 
(circles) and parabolic and lattice potential found analytically (stars). Shown are only 
the first few eigenvalues for f2 = 0.1 and Vo = 0.3. A similar result is demonstrated 
in panel (b) but for a larger lattice depth, Vq = 0.5 (and the same value of fi). For 
the latter case (Vq — 0.5), panels (c), (d) and (e) show the first few eigenmodes for the 
parabolic potential (thick solid), parabolic and lattice potential computed numerically 
(thin solid), and the same ones given by Eq. (fl4|) (dashed). The potential (rescaled for 
visibility) is shown by the dash-dotted line. 



using a transverse trap frequency u± = 2tt x 5Hz, the parameter ft is of order of 10 -4 ; 
thus, it is a natural small parameter of the problem. 

Our analysis starts by considering the non- interacting limit s — > 0, in which the 
GP equation becomes a linear Schrodinger equation. Then, seeking stationary localized 
solutions of the form ip(x,y,t) = exp(—iE mtn t)u mjn (x,y) [where E m ^ n are discrete 
energies and u mtn (x,y) are the corresponding linear eigenmodes], and rescaling spatial 
variables by y/Ti, we obtain the following equation: 



E„ 



-u 



m,n 



1 . 

-Am. 
2 



,,,.„ + -i-'- 2 • !f )n, 



+ 



Vo 
ft 



cos 




(6) 



+ cos 



2y 



it, 



The next step is to separate variables through u(x, y) = u m (x)u n (y) and split the energy 
into E m>n = E m + E n , to obtain two ID eigenvalue problems of the same type as above, 
namely, 



1 d u m 1 n 

2^ + r u ' 




E n 
"ft 



■u, 



(7) 



and a similar one for y (with x replaced by y and the subscript m replaced by n) 
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We now restrict ourselves to the physically relevant regime of < <C 1 as 
discussed above. In this case, we may use v = Vfi as a formal small parameter 
and develop methods of multiple scales and homogenization techniques [16] in order 
to obtain analytical predictions for the linear spectrum. In particular, introducing the 
fast (i.e., rapidly varying) variable X — xjv and rescaling the energy as e m = E m /Q, 
the eigenvalue problem of Eq. (J7|) is expressed as follows: 

U 2 C U ~ V 7TT^? + ^OL ) Urn = ^Vmi (8) 



where 



dxdX 
1 d 2 

i d 2 



1 X 

2 Ox 2 2 



OL 



2<9X 2 



+ y cos(2X) 



(note that Vq/VI was treated as an 0(1) parameter) 
series expansion (in v) for u m and e m [16j, namely, 



^ 71 1 



U + VUx + V U2 + ... 

E- 2 V~ 2 + E-iV~ l + e 



E\V + . 



(9) 
(10) 

Additionally, we consider a formal 



(11) 
(12) 



To this end, substitution of this expansion in the eigenvalue problem of Eq. (JHJ) and 
use of the solvability conditions for the first three orders of the expansion [i.e., 0(1), 
0{y) and 0{v 2 )] yields the following results for the eigenvalue problem of the original 
operator. The energy of the m-th mode can be approximated by: 



En 



1 



1, 



4 



(13) 



while the corresponding eigenfunction is given by: 

x 



U m [X) 



X 




(14) 



where a 

J2 



= (2 rn m!-\/7r)~( 1 / 2 ' ) is the normalization factor and H m (x) = 
e" {— l) m (d m / dx m )e~ x2 are the Hermite polynomials. 

Combining the results in the two orthogonal directions x and y yields a total energy 
eigenvalue 



E = --Ur \ (1- 2 V o) ^ - 1 



(15) 



and a corresponding eigenfunction which up to normalization factors can be written as: 



u r . 



(x, y) oc H, 



x H, 
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x exp 



( 



2 



2 



) 



(16) 



The particularly appealing feature of this expression is that it allows us to combine 
various (ground or excited) states in the x direction with different ones along the 
y direction. In this work we will focus on the simplest possible combinations of 
m, n G {0, 1, 2} and examine the various states generated by the combinations of these, 
which we hereafter denote \m,n). In particular, below we will focus on the ground 
state 1 , 0) and the excited states |1,0), |1,0) + |0, 1), |1, 1), and |2, 0). Our aim is 
to investigate which of these states persist in the nonlinear regime in both cases of 
attractive and repulsive interatomic interactions, and study the stability of these states 
in detail. Notice that we will only illustrate (by an appropriate curve in the numerical 
results that follow) the linear limit, E mi7l of Eq. ( IT51) above, for the case of attractive 
interactions. 

3. Numerical Results 

3.1. The non-interacting limit 

We start by examining the validity of the above analytical predictions concerning the 
linear limit of the problem, namely, Eqs. ( |T5l) and ( jTBT) . The results are summarized 
in Fig. [H Panel (a) shows the ID harmonic oscillator energy spectrum (for Q = 0.1) 
and compares it with the energy spectrum obtained from numerical and approximate 
theoretical [see Eq. ffl3]) ] solutions of the combined harmonic and optical lattice potential 
for Vq = 0.3. Panel (b) offers a similar comparison but for a larger lattice depth, 
namely Vq = 0.5. One can clearly see that the theoretical calculation approximates very 
accurately the numerical results for the first few states (i.e., n — 0, 1, 2), while deviations 
become more significant for higher-order excited states. Panels (c), (d) and (e) show the 
zeroth, first and second eigenfunction of the purely harmonic potential (thick solid line), 
as well as of the harmonic trap and optical lattice potential as found numerically (thin 
solid line) and analytically [given by Eq. (JHJ)]. The green dash-dotted line represents 
the form of the combined potential. We once again note the good agreement of our 
analytical results above in comparison with the full numerical computation. 

3.2. The approach for the interacting case 

We now consider the full nonlinear problem of Eq. ([3]). In the following, we will monitor 
the two-dimensional, Vq x /x, parameter space (where fi is the chemical potential of the 
relevant modes) for nonlinear excitations that stem from the linear spectrum of the 
problem. We perform the relevant analysis first in the case of attractive interatomic 
interactions and then in the case of repulsive ones. Notice that the parameter Q will 
be fixed to a relatively large value, namely Q = 0.1; this is done for convenience in 
our numerical simulations (such "large" values of Q correspond to smaller condensates 
that can be analyzed numerically with relatively coarser spatial grids), but we have 
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checked that our results remain qualitatively similar for smaller values of Q (results not 
shown here). Nevertheless, it should be noted that even such a value of Q, together with 
the considered range of values of the other normalized parameters (chemical potential, 
number of atoms, lattice depth, etc - see below) is still physically relevant. For example, 
our choice may realistically correspond to a 87 Rb condensate containing ~ 15, 000 atoms, 
confined by a harmonic potential with frequencies uj z = 20uj± = 2tt x 240Hz (so that 
uj± = 2tt x 12 Hz) and an optical lattice potential with a periodicity d ~ 3/im. The 
recoil energy in this case is E rec /h = 60 Hz (so that ujl = 2ti x 120 Hz, giving Q = 0.1), 
and a lattice depth of V = 0.3 corresponds to ~1.2 E TCC . To further set the scale for 
the simulations described below, such a BEC with repulsive interactions in the purely 
harmonic trap (where the optical lattice is not applied) would have a chemical potential 
of /i ~ 0.5 in our dimensionless units. 

In the following sections, stationary solutions of the full nonlinear problem are 
sought in the form 

ip(x, y, t) = exp(-z/i min t)w min (x, y), (17) 

where ^ m , n (which is the nonlinear analog of the energy E m ^ n found in the non-interacting 
limit) represents the chemical potential. Note that we will henceforth avoid using 
subscripts m and n when the meaning is clear, in the interest of avoiding notational 
clutter. 

3.3. Attractive interatomic interactions 

3.3.1. Existence and Stability We begin by looking at the case of attractive interatomic 
interactions. The most fundamental solutions are those belonging to the |0, 0) branch, 
which represents the ground state of the system and is shown in Fig. [21 The top left panel 
of this figure shows the diagnostic that we will typically use to follow the V x /i surface, 
namely the rescaled number of particles N(Vo,fi) = J \u m ^ n \ 2 dxdy as a function of the 
chemical potential /i introduced above, and the optical lattice depth Vq. Essentially, 
the grayscale values in this plot correspond to the number of atoms needed to obtain a 
particular chemical potential with a particular lattice depth; lighter values correspond 
to more atoms. As N becomes smaller, through the appropriate variation of /i, we 
approach the linear limit so one expects the solution to degenerate to the corresponding 
linear eigenmode (for [i tending to the corresponding eigenvalue of the linear problem). 
Figure [2] shows our observations for this fundamental branch, which seems to disappear 
for /i mi „(V r ) « E m , n (V ) = —\Vq + (l — \V{f*j Vt (n + m + 1). Naturally the surface 
degenerates to its linear limit for /i m) „(Vo) — > E mtn (y ) and the number of particles is a 
decreasing function of /i, contrary to what is the case in the repulsive nonlinearity (see 
below). This is a well-known difference between the two cases that has been documented 
elsewhere (see, e.g., Refs. [HI [16]). It is relevant to indicate here that although there is 
a direct correspondence between the atom number N and the chemical potential /i, we 
opt to illustrate our results as a function of /x (and Vo), since the latter is the relevant 
parameter entering the mathematical setup of the problem and it is the one for which 
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we developed an analytical prediction in the linear limit. Nevertheless, we also give 
N(Vo,fi), so as to associate in each case the relevant chemical potential (and lattice 
strength) for a given configuration with the corresponding physical quantity, i.e., the 
atom number. 

It is important to highlight here that that the numerical computations have been 
performed in a domain of size 201 x 201, with Ax = Ay = 0.15. The size of the 
grid weakly affects the value of the respective eigenvalues. In particular, the energy 
eigenvalues in the non-interacting limit of the one-dimensional decoupled eigenvalue 
problem (with Vq — 0.5) for the n = and n = 1 mode that we report below are 
Eq m —0.01453 and E\ w 0.07639 respectively for this coarser domain, while for the a 
domain of size 3001 with Ax = 0.01 they become —0.01408 and 0.07692 respectively 
(computations show the discrepancy is uniformly smaller for smaller values of Vo). This 
feature will weakly affect the quantitative aspects of the results that follow (in essence, 
providing an error bar in the estimates below of A/i « 5 x 10~ 4 and similar for the 
eigenvalues A introduced below), but is essentially necessary, given the limitations of 
standard eigenvalue solvers for such big Jacobian eigenvalue problems. 

The linear stability of the solutions is analyzed by using the following standard 
ansatz for the perturbation, 



where A = A r + zA, is the generally complex eigenvalue (subscripts r and i denote, 
respectively, the real and imaginary parts of A) of the ensuing Bogoliubov-de Gennes 
equations [TJ [2} [3l H] , and (a, b) T is the corresponding eigenvector. The real part A r of the 
eigenvalue then determines the growth rate of the potential instability of the solution, 
since for positive real values the perturbation will grow exponentially in time. The 
imaginary part \ denotes the oscillatory part (frequency) of the relevant eigenmode. 
The top right of Fig. [2] depicts the stability domain, denoted by S{V ,fi) = max A (A r ), 
in terms of the maximum real part of all eigenvalues as a function of the lattice depth 
Vo and the chemical potential \i. This quantity S is a particularly important one from a 
physical point of view since if a perturbation to the system has initially a projection p(0) 
onto the most unstable eigenmode of the linearization, then this perturbation will grow 
in time according to = | \p{0) 1 1 exp(St) while the solution is sufficiently close in 

space to the original profile. Hence, given the initial condition profile [which determines 
p(0)] and S, we can determine for an unstable configuration the time t required for 
the instability to manifest itself, i.e., for p{t) to become of the order of the solution 
amplitude. 

It is important to note, in connection to our numerical linear stability results, that 
the |0, 0) branch can become unstable for // < /i cr (Vo) (see top right panel of Fig. [2]) due 
to the appearance of a real pair of eigenvalues. This instability for large N is something 
that may be expected in the case of attractive interactions under consideration, as the 
corresponding 2D GP equation for an a homogeneous BEC (i.e., without any external 
potential) is well-known to be subject to collapse [25J. However, it should also be 




= e -W [ u (x, y) + (a(x, y)e xt + b*(x, y)e yt )\ , 






Figure 2. (Color Online) The ground state in the case of attractive interatomic 
interactions. The top left panel shows the rescaled number of particles N(Vo,fi) = 
J \u\ 2 dxdy as a function of the amplitude of the optical lattice Vq and the chemical 
potential [jl; the red line represents the approximation of the energy eigenvalue E(Vo) 
of the linear problem given by Eq. (fT5|) . For each Vq, the number of atoms, Nv (/J.) 
approaches zero in the limit fi — *■ E(Vq). The top right panel shows the stability domain 
S(Vo, ix) = max(A r ); the red line here corresponds to the stability window S < 10~ 4 . It 
is clear that for each Vb, there is a window of values of /i for which Sv (/i) < 10~ 4 . This 
is expected, since the attractive nature of the interatomic interactions leads to blowup 
above a critical value of N. The bottom left and right panels depict, respectively, 
a contour plot of an unstable solution u in the [x, y) plane and its corresponding 
spectral plane (A r ,Aj) [for (Vo,fi) = (0.21,-0.23) corresponding to the parameter 
value depicted by the red circle in the panels of the top row] . The bottom-left colorbar 
represents atomic density. In the bottom-right plot, the presence of real eigenvalue 
pairs denotes instability (its growth rate is given by the magnitude of the real part), 
while the imaginary pairs indicate the frequencies of oscillatory modes. 
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Figure 3. (Color Online) Similar to Fig. [2] but for the case of the |1,0) state for 
attractive interactions. The results shown in the bottom row correspond to parameter 
values (Vb,£t) = (0.21, —0.081) (see red circle in top panels). 



expected that very close to the linear limit the growth rate of the instability is essentially 
zero (cf. with the top left panel of Fig. 2 of Ref. [TTJ for Vo = which is not shown here). 
Essentially, the potential appears to stabilize the solitary wave against dispersion in this 
regime (i.e., close to the linear limit), but cannot stabilize it against the catastrophic 
collapse-type instability. Furthermore, in the presence of the optical lattice we can 
observe that there is always a range of chemical potentials for which the condensate is 
stabilized, in accordance with what was originally suggested in Ref. [26]. Furthermore, 
even in the 3D case it is in principle possible to arrest collapse by appropriate choices 
of the parameters [27] . 

Next, we consider real- valued solutions with m + n — 1. The |1,0) state (again in 
the case of attractive interactions) is shown in Fig. [3l This branch is always unstable, 
due to up to two real eigenvalue pairs and one complex quartet. A typical example of 
the branch in the bottom panels of the figure reveals this instability. 

The 1 1,0) + |0, 1) configuration for the attractive interactions case is shown in 
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Figure 4. (Color Online) The state |1, 0) + |0, 1) for attractive interatomic interactions. 
The top row is similar to Fig. [TJ the middle row is for parameter values (Vb,/i) = 
(0.31, —0.201) (see blue cross in top panels), and the bottom one is for (0.35, —0.481) 
(see red circle in top panels). 
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Figure 5. (Color Online) The state |1, 1) for attractive interatomic interactions. The 
layout of the figure is similar to the one used in the previous figures. The parameters 
for the solution depicted in the middle and bottom rows are (Vb, /i) = (0.2, —0.181) (see 
blue cross in top panels) and (0.45,0.039) respectively (see red circle in top panels). 
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Figure 6. (Color Online) The state |2, 0) for attractive interatomic interactions. The 
layout of the figure is similar to the one used in the previous figures. The parameters 
for the solution depicted in the bottom row are (Vq, fi) = (0.3, —0.001) (see red circle 
in top panels). 



Fig. HI This configuration turns out to be unstable in a large fraction of the regime of 
parameters considered due to a quartet of complex eigenvalues. However, remarkably, 
as Vq and \i are increased and decreased respectively, it is possible to actually trap this 
state in a linearly stable form (eliminating the relevant oscillatory instability). This 
indicates that it would be of particular interest to try to identify such a state (which 
resembles an out-of-phase soliton pair) in a real experiment. Also, as expected, the 
solution degenerates to its linear counterpart as \i — > E. Images of a typical unstable 
solution and its complex quartet are shown along with a stable solution from the top 
right-hand region of the two-parameter diagram. 

We should also note in passing that states in the form of |1, 0)+i|0, 1) would produce 
a vortex waveform; however since such states have been studied in some detail earlier 
in Ref. [21] in a similar setting (i.e., in the presence of an external potential containing 
both harmonic and lattice components), we do not examine them in more detail here. 
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(a) (b) 




(c) (d) (e) 




Figure 7. (Color Online) Dynamics of the unstable states (in the case of attractive 
interatomic interactions) that were shown in the previous figures. Shown are space-time 
evolution plots given by a characteristic density isosurface Dk = {x, y, t \ \u(x, y, t)\ 2 = 
k}, where k = a(max{ x y }{\u(x, y, 0| 2 }) and [a=0.7 for (a), 0.5 for (b) and (d), and 
0.3 for (c) and (e)]. (a) Ground state 1 0, 0) , which collapses very quickly, (b) First 
excited state |1,0), which collapses shortly after the ground state, (c) Degeneration 
of a 1 1,0) + |0, 1) state into an eventual single-pulse structure that survives for a long 
time after the merger. The unstable |1,1) (d) and |2, 0) (e) states deform, for very 
short times as expected from the strong instabilities identified in their spectra, and 
subsequently collapse. 

We now turn to solutions featuring m + n = 2. First, we consider the |1, 1) branch 
for the attractive case in Fig. [5j In this case the solution may possess between one and 
three complex eigenvalue quartets in its linearization (the middle panel of the figure 
shows a particular unstable case where there are two such quartets). However, once 
again, there exists a region in the right side of the relevant parameter space [i.e., for 
appropriate (Vo,/x)] where the solution is found to be linearly stable and all potential 
oscillatory instabilities are suppressed. The bottom panel of Fig. [5] shows such a linearly 
stable case of the quadrupolar configuration, which, again, should be experimentally 
accessible. 

Finally, we consider the state |2,0), as depicted in Fig. [6j This configuration is 
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highly unstable throughout our parameter space, with up to four real pairs and one 
complex quartet of eigenvalues. A typical example of the unstable configuration and its 
spectral plane of eigenvalues is shown in the bottom panel of Fig. EJ 

3.3.2. Dynamics Now we corroborate our existence and stability results (for the 
attractive interactions case) with an investigation of the actual dynamics of typical 
unstable solutions selected from the above families. For each case, the particular solution 
presented in the corresponding figures is perturbed with random noise distributed 
uniformly between —0.05 and 0.05 and integrated over time. It is important to note 
that although a random perturbation is used here to "emulate" the experimental noise, 
the system is deterministic and the sole relevant feature of any (generic) random 
perturbation is its projection onto the most unstable eigendirection(s) of the perturbed 
solution profile. These projections, as indicated above, will grow (determinstically) 
according to the corresponding growth rate. For the time propagation, we implement a 
standard 4th-order Runge-Kutta integrator scheme where we have numerical consistency 
and stability for the conservative time step of At = 10 -3 . The results are compiled in 
Fig. El 

Panel (a) in Fig. [7] depicts the catastrophic instability of the ground state, |0,0), 
which is subject to collapse, occurring at t ~ 25. Panel (b) shows similar behavior for 
the |1, 0) state, in which the two lobes appear to self-focus independently, although one 
eventually prevails and collapses for t « 35. It is very interesting to note that while the 
|1, 0) state collapses, its more stable superposition with the |0, 1) state survives for longer 
times, as expected, and also eventually merges into a ground-state-like (single pulse) 
configuration, which was found to have a number of atoms just on the unstable side 
of the boundary of stability for such a structure. The resulting state actually survives 
for very long times, oscillating within one of the wells where it originally collected itself 
[see the panel (c)], apparently stabilized by the ensuing oscillations. Panels (d) and (e) 
show, respectively, the relatively rapid break up and subsequent collapse of the |1,1) 
and |2, 0) states. 

3.4- Repulsive interatomic interactions 

3.4-1. Existence and Stability Now we will investigate the results pertaining to 
repulsive interatomic interactions for the same linear states examined above. We once 
again start with the ground state |0,0) branch shown in Fig. [8l The top left panel of 
Fig. M shows the number of atoms iV(Vo,/i) as in the previous section. However, the 
linear stability S(V , /i) for this case is omitted because, as may be expected, this branch 
is stable throughout the parameter space, in contrast to its attractive counterpart (which 
is subject to collapse). 

We now turn to excited states with m + n = 1. Figure [9] shows features similar to 
the previous one, but now for the state |1, 0). This branch is found to always be unstable 
due to the appearance of up to three real eigenvalue pairs. The top panels depict the 
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Figure 8. (Color Online) The ground state state |0, 0) solution in the case of repulsive 
interatomic interactions. The bottom panels correspond to (Vo,fj) — (0.3,0.37). The 
stability of the ground state persists over the entire parameter space, and hence the 
stability surface is omitted from this set. 



dependence of the number of atoms N(V , fi) (left), and instability growth rate S(V , //) 
(right) on the lattice depth V and the chemical potential \i. A sample profile and its 
eigenvalue spectrum are given in the bottom panels, indicating the presence, in this 
case, of three real eigenvalues pairs. 

The next state we consider is the |1, 0) + |0, 1) state which is presented in Fig. [TUJ 
This state always possesses a quartet of complex eigenvalues, and up to two additional 
pairs of real eigenvalues, and is thus unstable for all fi. It it worth noticing, however, 
that the instability weakens to relatively benign small magnitude complex quartets for 
intermediate values of the chemical potential, roughly /i e (0.4,0.9), and large lattice 
depths, Vo > 0.3. This suggests that such a configuration should be long-lived enough 
that it could be observable in experiments with repulsive condensates. 

Next, we consider the states with n + m = 2 (again for repulsive interatomic 
interactions), starting with the |1, 1) branch in Fig. [TTJ The branch is also always 
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Figure 9. (Color Online) The state |0, 1) in the case of repulsive interatomic 
interactions. The bottom row illustrates a sample profile (left) and the corresponding 
eigenvalue spectrum (right) for (Vo,u) — (0.2,0.47) displaying the strong instability 
arising from the three real eigenvalue pairs. 

unstable, possessing a complex quartet and then up to four additional real pairs for 
larger values of fi. The instability is shown in the right subplots of Fig. [Ill where 
the spectral plane of the bottom right panel corresponds to the solution of the bottom 
left one, for parameter values (Vo,/x) = (0.2,0.61). It is interesting to note that such 
states are reminiscent of the domain-walls presented in Ref. [28j (here the domain-wall 
is imposed by the difference in phase), which however were found as potentially stable 
structures in multi-component condensates. 

Next, the case of the |2, 0) state is shown in Fig. [12j Here, there are up to three 
complex quartets along with three real pairs of eigenvalues. It is notable that for higher 
values of the lattice depth, these states are deformed as the lattice "squeezes" the central 
maximum separating the two minima (see middle left panel of the figure). The contour 
plot shown in the bottom right panel suggests that further increase of lattice depth may 
lead to a new configuration altogether when the two local maxima eventually pinch off of 



Structure and stability of 2D BECs under harmonic and lattice confinement 



20 



0.21 
0.4 
3-0.6 
0.8 
1 



0.1 0.2 



0.3 









0.5 






0.6 






0.4 






0.2 






^ 






-0.2 






-0.4 




■ 


-0.6 

-0.5 



O § 



o § o 



0.2 
0.15 
0.1 
0.05 



-0.05 



, 0.05 

A 



0.1 




0.5 



M~ 



-0.5 











o 
o 






° § ° 


o 
o 









o 
o 


o 1 

( 
1 


i 

) 
) 


o 













-0.01 



, 0.01 

K 

r 



Figure 10. (Color Online) Same as the previous figures but for the |1, 0) + |0, 1) state 
in the case of repulsive interatomic interactions. This state is always unstable due to 
at least an eigenvalue quartet and up to two other real pairs. Note that there exists 
a region of weak instability, where only small magnitude quartets are present. The 
middle and bottom panels show the contour plots of this state and its linearization 
spectrum for (Vb,/i) = (0.2,0.55) (see blue cross in top panels) and (0.37,0.59) (see 
red circle in top panels), respectively. 
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Figure 11. Same as in Fig. [SI but for the state |1, 1) in the case of repulsive interatomic 
interactions for parameter values (Vb,/x) = (0.2,0.61). 

each other. This deformation is a direct consequence of the presence of the (repulsive) 
nonlinearity, which results in drastically different configurations as compared to the 
linear limit of the structure. 

3.4-2. Dynamics We performed numerical simulations to investigate the evolution of 
typical unstable states in the case of repulsive interatomic interactions, using similar 
time-stepping schemes as discussed above in the case of attractive interactions. Apart 
from the ground state, all excited states presented in the previous section were predicted 
to be unstable, and this is confirmed in this section. In the particular case of the state 
|1,0) + |0, 1) which was found to be weakly unstable (see bottom row of Fig. [10]), the 
instability takes a considerable time to manifest itself. The evolution of this state is 
depicted in panel (b) of Fig. [131 it is clearly seen that the state persists up to £ ~ 300. 
On the other hand, panel (a) shows the evolution of the state |1, 0) which persists up to 
t ~ 40, while the bottom row of panels shows the dynamics of the (c) |1, 1) state and of 
the (d) 1 2 , 0) , both persisting also up to t « 40. All of these excited states degenerate 



Structure and stability of 2D BECs under harmonic and lattice confinement 



22 



0.3 
0.4 
0.5 
0.6 
0.7 






15C 


0.3I 






0.4 












0.5 










50 


0.6 








0.7 




0.25 

0.2 

0.15 

0.1 

0.05 



0.1 0.2 0.3 0.4 

V„ 



0.1 0.2 0.3 0.4 

V„ 




15 
10 
5 

-5 
-10 
-15 



-10 



-10 





x 





x 



10 



-0.1 -0.05 . 0.05 0.1 

A 




10 





15 


0.4 


10 


0.2 


5 





>• 


-0, 




-0.4 


-5 


-0.6 


-10 
-15 




-10 





x 



10 



0.5 



|~0.5 



Figure 12. (Color Online) Same as in Fig. [9) but for the state |2, 0) (in the case of 
repulsive interatomic interactions) with parameter values (Vo,/i) = (0.2,0.7) (see red 
circle in top panels). The bottom row shows profiles for smaller, Vq — 0.1 (left, see blue 
cross in top panels), and larger, Vq — 0.35 (right, see green diamond in top panels), 
values of the optical lattice depth. 
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Figure 13. (Color Online) The dynamics of the unstable states in the case of repulsive 
interatomic interactions. Panels (a) and (b) show, respectively, the evolution of the 
states |1,0) and |1,0) + |0, 1). It is clear that the state |1,0) + |0,1) is subject to a 
weaker oscillatory instability for the parameter values mentioned in the bottom row of 
Fig. 1101 and, as a result, the original configuration persists for a long time. Panels (c) 
and (d) show the dynamics of the states with m + n = 2, namely (c) |1, 1), and (d) 
1 2, 0) . All these solutions ultimately degenerate into ground-state-like configurations. 
The density isosurfaces are taken at a = 0.5 with the exception of (d) at a = 0.4. 

into ground-state-like configurations. Notice that transient vortex-like structures seem 
to appear during this process but they do not persist in the eventual dynamics and are 
hence not further discussed here. 

4. Conclusions and discussion 

In summary, we have studied the structure and the stability of a pancake-shaped 
condensate (with either attractive or repulsive interatomic interactions) confined in a 
potential with both a harmonic and an optical lattice component. Starting from the 
non- interacting limit, and exploiting the smallness of the harmonic trap strength, we 
have employed a multiscale perturbation method to find the discrete energies and the 
corresponding eigenmodes of the pertinent 2D linear Schrodinger equation. Then, we 
used the results found in this linear (non-interacting) limit in order to identify states 
persisting in the nonlinear (interacting) regime as well. This investigation revealed that 
the most fundamental states (emanating from combinations of the ground state and 
the first few excited states in the two orthogonal directions of the optical lattice) can 
indeed be continued in the nonlinear regime. To demonstrate this continuation, we used 
two-parameter diagrams involving the effective strength of the nonlinearity (through 
the chemical potential) and the optical lattice depth. 

Excited states were typically found to be unstable. The instability was found to 
result in either wavefunction collapse or a robust single-lobed structure in the case of 
attractive interactions; on the other hand, in the case of repulsive interactions, the 
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instability was always found to lead to the ground state of the system. Nevertheless, 
noteworthy exceptions of stable or very weakly unstable states were also revealed. These 
include the |1,0) + |0, 1) and the |1,1) states in the case of attractive interatomic 
interactions. Moreover, in the case of repulsive interactions, the same state, |1, 0) + |0, 1), 
was found (in certain parameter regimes) to be only very weakly unstable. Direct 
numerical simulations confirmed that the instability of this state is indeed weak, and it 
manifests itself at large times, an order of magnitude larger than the ones pertaining 
to the manifestation of instabilities of other excited states. Thus, it is clear that the 
obtained results suggest that the state |1, 0) + |0, 1) has a good chance to be observed in 
experiments with either an attractive or a repulsive pancake condensate. It is especially 
important to highlight that these states are stabilized (or quasi-stabilized) only in the 
presence of a sufficiently strong optical lattice; hence the latter potential plays a critical 
role in determining the stability of the states presented herein. 

As described in Section III B, the parameters used in our analysis have been 
chosen in order to facilitate convenience of the numerical computations, while also 
within range of experimentally achievable limits of atom number, chemical potential, 
harmonic oscillator frequencies, and optical lattice depth and periodicity. Furthermore, 
we believe that the stability and spatial structure of the states examined here can be 
examined experimentally. For example, we imagine utilizing a BEC held in a pancake- 
shaped harmonic trap, created by an optical field. By using an optical trap rather than 
a magnetic trap, the scattering length of a BEC may be adjusted using a Feshbach 
resonance. We envision that an optical lattice potential is ramped on and superimposed 
on a BEC with an interatomic scattering length tuned to be near zero. Once the lattice 
has reached the desired depth, the scattering length can be further adjusted with a 
magnetic field to be either positive or negative (the latter option would need to be 
within a region of stability that does not result in collapse of the BEC). Finally, phase 
imprinting techniques [29] can be used to generate the desired phase profile of the BEC. 
By optically examining the state of the BEC at various points in time after phase profile 
imprinting, the stability of the generated states can then be examined and compared 
with our numerical results and stability analysis. For example, with an optical lattice 
frequency of ujl = 2ir x 120 Hz (as in the example of Section III B), the time unit of our 
dynamical evolution plots is 1.3 ms. This implies that for the cases we have examined, 
signatures of instability would be typically visible on the experimentally feasible 10 to 
100 ms timescale. We therefore believe that our predictions could be examined with 
current experimental techniques. 

There are various directions along which one can extend the present considera- 
tions. A natural one is to extend the analysis to fully 3D condensates and examine 
the persistence and stability of higher- dimensional variants of the presented states. A 
perhaps more subtle direction is to consider a different basis of linear eigenfunctions 
in the 2D problem, namely instead of the Hermite-Gauss basis used here, to focus on 
the Laguerre- Gauss basis of the underlying linear problem with the parabolic potential. 
Under such a choice, it would be interesting to examine how solutions of that type, in- 
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eluding one-node and multi-node ring-like structures (see, e.g., Ref. [T7j and references 
therein), are deformed in the presence of the lattice and how their stability is corre- 
spondingly affected. Finally, as discussed above, it appears that the setting considered 
herein should be directly accessible to present experiments with pancake-shaped BECs. 
In view of that, it would be particularly relevant to examine which ones among the 
structures presented in this work can survive for evolution times that are of interest 
within the time scales of an experiment. 
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